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ABSTRACT 

One of the most interesting and challenging aspects of formation guidance law design is the coupling of the orbit 
design and the science return. The analyst’s role is more complicated than simply to design the formation geometry 
and evolution. He or she is also involved in designing a significant portion of the science instrument itself. The 
effectiveness of the formation as a science instrument is intimately coupled with the relative geometry and evolution of 
the collection of spacecraft. Therefore, the science return can be maximized by optimizing the orbit design according 
to a performance metric relevant to the science mission goals. In this work, we present a simple method for optimal 
formation guidance that is applicable to missions whose performance metric, requirements, and constraints can be 
cast as functions that are explicitly dependent upon the orbit states and spacecraft relative positions and velocities. 
We present a general form for the cost and constraint functions, and derive their semi-analytic gradients with respect 
to the formation initial conditions. The gradients are broken down into two types. The first type are gradients of 
the mission specific performance metric with respect to formation geometry. The second type are derivatives of the 
formation geometry with respect to the orbit initial conditions. The fact that these two types of derivatives appear 
separately allows us to derive and implement a general framework that requires minimal modification to be applied to 
different missions or mission phases. To illustrate the applicability of the approach, we conclude with applications to 
two missions: the Magnetospheric Multiscale mission (MMS), and the Laser Interferometer Space Antenna (LISA). 

INTRODUCTION 

Developing guidance law algorithms for spacecraft formations and constellations is a challenging problem that can 
seemingly involve mutually exclusive goals. From the perspective of a mission analyst working on orbit design for 
distributed spacecraft missions, it is desirable to have a flexible method that can be applied to many problems, and 
that permits modification of the cost and constraints with as little effort as possible. We also desire an approach that 
can handle real-world mission constraints, find an optimal solution in a reasonable amount of time, and converge 
consistently with relatively poor initial guesses. The work presented in this paper was motivated by the need for 
an algorithm that can provide formation initial conditions for many different phases in the mission design process. 
These phases include, but are not limited to: developing initial conditions to use as target states to achieve after 
separation from the launch vehicle and separation from the spacecraft stack, developing new conditions for formation 
maintenance after perturbations and errors in navigation and control have caused a significant degradation in mission 
return, and developing new formation initial conditions for a new mission phase when the desired formation geometry 
is significantly different. For each of these mission phases we need to determine initial conditions for the orbits of 
multiple, cooperating spacecraft, that optimally satisfy mission requirements and constraints. This is a long list of 
goals, but we can come very close to meeting them all by developing effective design algorithms. 

The approach presented in this paper is a direct parameter optimization method which assumes that at a given 
instant in time, one can formulate a performance metric that is an explicit function of the inertial cartesian states, 
and the relative position vectors, relative velocity vectors, range and range rates of all spacecraft in formation. The 
performance metric allows us to quantify the effectiveness of the relative dynamics of the spacecraft at a particular 
point in the orbit. The cost function, J, is simply the integral of the instantaneous performance metric over regions 
of interest to the particular mission. The goal is to maximize or minimize the cost function, according to the specific 
mission, while simultaneously satisfying equality and inequality constraints consistent with mission requirements and 
real-world limitations. 

We begin by Dosing a continuous-time set of cost and constraint functions that define a general framework for the 
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problem. The continuous-time system is the problem we would like to solve. However, there are certain advantages 
to discretizing the problem. Hence, we discretize the problem to enable a Nonlinear Programming (NLP) routine 
to solve for an optimal solution. The discrete cost function can be calculated by quadrature, or a simple trapezoid 
rule, using numerical integration of the non-linear orbit equations of motion, with perturbations, to give us the 
state of the formation at discrete times in the orbit. The constraints are calculated in a similar manner. We show 
how to calculate the derivative of the cost function and the Jacobian of the constraint functions using the State 
Transition Matrix (STM). The STM is found by numerically integrating an additional 36 equations per spacecraft 
with terms that include perturbations. There is a significant advantage to this approach. Numerical integration, 
while non-trivial, is a solved problem and something that modern computers perform efficiently. Secondly, as we 
will see in a later section, the derivatives of the cost and constraints contain two types of terms. The first type of 
terms are composed of portions of the STM and inertial states. The second type of terms contain derivatives of 
algebraic functions with respect to the formation geometry. Because the terms are separated in this way, we can 
implement the algorithm in a very general way so as to allow minimal modification to new problems. We conclude 
the paper with two applications that illustrate how the method is applied to real-world problems. The first example 
is the Magnetospheric Multi-Scale (MMS) mission. The second example is the Laser Interferometer Space Antenna 
(LISA). Let’s begin by defining some notation. 

SYMBOLS 

Variables 

r Position vector 

v Velocity vector 

n Number of spacecraft 

m Number of unique sides in the formation 

N Number of path constraints 

rifc Number of quadrature points 

Ck Quadrature constant at point k 

Sik Vector defining side z at quadrature point k 

s ik Rate vector of side z, at point k 

s t k Length of side z, at point k 

Rate of change side z, at point k 
Orbit state transition matrix 
A Upper left 3x3 partition of $ 

B Upper right 3x3 partition of 

C Lower left 3x3 partition of 

D Lower right 3x3 partition of 

Subscripts 

z Side index 

k Quadrature point index 

j Spacecraft index 

£ Spacecraft index 

p Dummy index 

J Indicates association with cost function 

C Indicates association with constraints 

General Problem Statement 

The reliability and effect iveness of an optimization algorithm is intimately dependent upon the parameterization of 
the cost and constraint functions. In this section, we pose a set of cost and constraint functions for a formation of 
spacecraft as explicit functions of the cartesian states of the spacecraft in formation, and the relative motion of the 
spacecraft. Let’s begin with some definitions. Recall that our goal is to find a set of initial cartesian states for a set 
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of n spacecraft. The vector defining the initial cartesian states for all spacecraft is 


X 0 = 


\>l 


T 

*4 


JT }T 


(i) 


where 

Xoi = [r^i vS] r = [Xoi Voi Zo 1 iol 2/ 0 i ioi] T (2) 

and so on. For our implementation, we choose to use the Earth MJ2000 equatorial coordinate system to express the 
cartesian states. However, it is possible to work in other coordinate systems, including rotating systems, as long as 
you are consistent. We assume that the equations of motion for the spacecraft are explicit functions of the spacecraft 
state and time, or: 

X(t) = f(X,t) (3) 

Let’s define the i th unique side vector of the formation, s*, as the vector from the £ th to the j th spacecraft 


Si(t) = Tj(t)-r e (t) 


( 4 ) 


The subscripts j and t can be chosen arbitrarily as long as for each value of i there is a given j, £ pair, the pairs 
result in a unique side, 1 < j < n and 1 < £ < n, and the chosen definition is used consistently throughout the 
implementation. Similarly to the side vector, the side rate vector s *, the side length s,, and the side rate Si are 
defined as 


S i{t) =Tj(t) -r e (t) 

( 5 ) 

<i{t) = (si(t) T Si(t)) 1/2 

(6) 

Si(t) = JJ) Si(t) 

( 7 ) 


Let’s pose a continuous-time cost function of the form 



Fj (X(t),Si(t), s i(t), Si(t), Si(t), Cj) dt 


( 8 ) 


where C j is a set of constants. Note that we assume the cost function is an explicit function of the relative geometry 
of the spacecraft. In Eq. (8), Fj is the performance metric that contains information on how well the geometry 
satisfies mission goals at a particular instant. Fj is mission specific and must be formulated according to a given 
problem. In a later section we present two examples for illustrative purposes. Similarly, we can define nonlinear 
equality constraints G^, and nonlinear inequality constraints, G /, as 


G e (X 0 ) = F E (X(t),Si(t),Si(t),Si(t) y Si(t),c E ) = C E (9) 

G/(X 0 ) = F / (X(t),s i (t),s i (t),s i (t),s i (t),c / ) < Cj (10) 

Finally, we also impose linear equality constraints, linear inequality constraints, and bound constraints, respectively 
as 

A*X 0 = B e (11) 

A/X 0 < B 7 (12) 

L < X 0 < U (13) 


We seek a solution that minimizes J, while simultaneously satisfying the constraints in Eqs. (9-13). However, 
there are several issues that make this difficult, if not impossible, to solve exactly. The first is that it is not possible 
to develop closed form solutions of X in terms of t. Similarly, we don’t know how to express the relative geometry it 
terms of explicit functions of time. Certain nonlinear constraint functions can also cause difficulties. By discretizing 
the problem, we can avoid these complications. Yet, we must be satisfied with a solution that is not exactly 
optimal. However, the difference between the rigorously optimal solution and the near-optimal solution obtained 
from parameter optimization is often so small as to make the difference academic. Now let’s take a look at how to 
discretize the continuous-time problem. 
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Discretization of the Cost and Constraints 


In general, we cannot solve the integral in Eq. (8). Also, it is often difficult to evaluate the constraint functions 
in Eqs. (9-13) along the entire trajectory. Furthermore, differentiating the equations with respect to X D is usually 
difficult. For these reasons, it is convenient to approximate the integral in the cost function by using a quadrature 
rule. Similarly, it is convenient to evaluate path constraints at discrete points along the trajectory. Once we have 
recast the cost and constraints as a discrete system, we can take the derivatives relatively easily, as we’ll see near the 
end of this section. 


Let’s begin by writing the cost function as a quadrature rule as follows: 

nk 

J (X 0 ) ~ Cj ^ ^ CfcJFj(Xfc, Sjfc, s j/-, Si ^ , Sik) (14) 

k = 1 

where k is the quadrature point index, rik is the number of quadrature points, Ck is the quadrature constant associ- 
ated with point fc, X^ = X(tfc), s ^ = s i(tfc) and so on. 


We have imposed many types of constraints in the continuous time system. Due to space limitations, we’ll limit 
the discussion here to the most difficult of these constraints: nonlinear path constraints. Nonlinear point constraints, 
linear, and bound constraints are much easier to handle, and we leave them as an exercise for the reader. The 
nonlinear path constraints are constraints that must be satisfied along a section of the trajectory and not at a 
specific point. We can group the equality and inequality path constraints together for convenience as 


G(Xo) 


G e (X 0 ) 

G/(X 0 ) 


( P\k (Xfc, Sjfc, Sj/s, Sik) ^ 
^2 k{X-k, S ik: S ifc , Sjfc, Sj/c) 

Ptk (X/j, Sjfc, S ik') s iki ik )) 


(15) 


\ ^iVfc(Xfe, S^fc, Sjfc, Siki $ik)) ) 


where the function Pek is the i th path constraint function evaluated at the k th quadrature point. Hence, we now 
have (N • k) point constraints in the discrete time system, whereas we only had N path constraints in the continuous 
time system. 


NLP algorithms are almost always more efficient when they are supplied with analytic derivatives as opposed to 
using finite differencing to find approximations for the derivatives. Recall that the parameters we allow the NLP 
algorithm to vary are the initial cartesian states, X G . We need to determine the derivative of the cost and constraints 
with respect to X 0 . We can write the derivative of the cost function as: 


8J 

0X o 


n k Q 

Cj / ^ Ck ~ v F /(Xfc, Sifc, S iki Siki $ik) 

ti ax ° 


Using the chain rule we can expand Eq.(16) to read 


dJ (dFjdXk dFjd s lk dFjds ik dFjds ik dFj 8s ik \ 

ox 0 ~ Cj 2^ Ck { dx k dXo + ds ik dXo ds ik ax 0 ds zk dx a + ds ik axj 

k—1 


(16) 


(17) 
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Similarly, for the constraints we can write 


dG 

dX 0 


( 


dP j k 

dX B 

dP2k 

ox 0 


(x fcl 

(x„, 


s ik j 


s ik ) 


dPtk 

0X o 


(X fc 




l 9Pnh 

\ dX 0 


(X fc 


s iki 


Sik: $ik) Si k ) 
^iki s iki ^i/e) 

Sj/cj Sik) Sik) 

Siki s ikj Sik)) 


\ 


( 18 ) 


Using the chain rule we can write the derivative of the I th constraint at the k th quadrature point as 

dPifc _ dPifc dX k dPtk ds ik 9P ek d$ik <9P^ fc dsik <9P ik ds ik 

dX Q ~ dX k 0X o ds ik 0X o ds ik dX Q ds^ dX 0 ds ik dX Q 


Inspecting Eqs. (17) and (19) we see that two types of terms appear and that these terms appear as pairs in a 
product. The first type of term is the derivative of the user supplied function with respect to the formation geometry. 
Examples of this type of constraint are 

dFj_ d^ dPtk 

0X k '' ds ik ’ ds ik { } 

The second type of term is the derivative of the geometry at a specific quadrature point, with respect to the formation 
initial conditions. Examples of this type of term are 


dX*k Os-ik &Sik /rtiy 

ax? ax? ax^ 1 ; 

What makes this method “simple” is two-fold. First, the terms of type one contain derivatives of the mission 
performance metric function with respect to the formation geometry. For many missions, the performance metric 
can be cast as an algebraic function of the formation geometry. As we’ll see in the examples presented in later sections, 
these derivatives are trivial. The second reason this method is “simple”, is that the second type of terms contain 
derivatives of the geometry with respect to the initial conditions, and the form of these derivatives is independent 
of the mission specific performance metric, and therefore they only have to be derived and implemented once! It 
is very convenient that the derivatives are separated in this way, and the fact that they are is simply due to the 
application of the chain rule to a cost function that is explicitly dependent upon the formation geometry. Let’s look 
at the derivatives of type two, which are the derivatives of the formation geometry and its evolution, with respect to 
the formation initial conditions. 


Derivatives of Formation Relative Geometry with Respect to Orbit Initial 
Conditions 


The derivatives of the formation geometry with respect to the formation initial conditions contain information on 
the sensitivity of the state of the formation at some time, with respect to the state of the formation at the 
initial epoch. We would expect portions of the STM to appear in these derivatives. Recall that the definitions for 
the geometry given in Eqs. (4-7) express the relative geometry, such as the relative position vector between two 
spacecraft, in terms of the cartesian states of the spacecraft. To calculate the derivative of side vector s with 
respect to the initial position of the p th spacecraft, r op , we can write 


ds^ 

Ot 0 p 


dr„ 


■ (r jk ~r ek ) = 


A j(tk)t 0 ) p = j 
-A e{t k ,t 0 ) p = £ 
0 otherwise 


( 22 ) 


where A j is the upper left 3x3 component of the STM for the j th spacecraft and so on. We define the four 3x3 
subcomponents of the STM as 

f \ — ( A3x3 (i) to) ^3x3 (tyt 0 ) ^ (23) 

V ^3x3(ht 0 ) 1^3x3 {hio) J 
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The derivative of s ik with respect to v op can be written as 

ds ik d 


dv op dv op 


(r jk - J 7fc) - 


B j(t k ,t 0 ) p = j 
-B e(tk,t 0 ) p — t 
0 otherwise 


The remaining derivatives of the formation relative geometry with respect to the initial conditions are: 

ds ik d 


dr op 

diiik 

dv op 


dr, 


(Vjk - v ek ) = 


op 


Cj(t k ,t 0 ) p = j 
C e (t k ,t 0 ) p = £ 
0 otherwise 


= ^r {Vjfc ~ Vtt) = 



_ 

d 

(s T s- t ) 1/2 

S T 

_ zih. 


dr op 

d r op 

\ s ik s ik) 

Sik 


ds^ 

d 

( s T uS ..) 1/2 

s T 
_ ^ik 


dv op 

dv 0 p 


Sik 

OSik 

1 

( -T 

T * s lk \ 

ds ik 

dr op 

Sik 

[■ s ik - 

S ik S ik 2 ) 

S ikS 

d r op 

ds ik 

1 

( -T 

T • S Ik 

ds ik 

dv„ p 

Sik 


S ik S ik o J 
S ik J 

dv op 


D j(t k ,t 0 ) p = j 
-T> t (tk,t 0 ) p = e 
0 otherwise 

d s ik 


vi op 

ds ik 


s Tk d*ik 
s ik dr 


op 


+ 


_ ik 


ds ^ 


Sik dv, 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


op 


The derivatives contained in Eqs. (22-30) only contain terms that involve the spacecraft states and their associated 
STMs. Therefore, these equations do not change from one problem to the next. What will change is the actual 
numbers in the STMs, and the spacecraft states. This enables us to implement the algorithm in a general way in 
order to minimize the amount of recoding to apply the algorithm to a new mission or mission phase. Now let’s take 
a look at some aspects of the implementation we use. Then, we’ll look at two mission applications. 


Implementation 

The implementation of a numerical optimization approach must be carefully considered to ensure speed and ease 
of use. There are several aspects of the implementation that will have an influence on the speed of the approach 
we present in this paper. The first is how the orbit states and STMs are propagated. It is recommended that all 
propagation is performed in a compiled language as opposed to an interpreted language such as MATLAB. This is 
primarily due to the fact that we must propagate 42 differential equations per spacecraft. For our implementation, 
we coded most of the algorithm in MATLAB, and use MATLAB ’s fmincon SQP routine. However, all propagation 
is performed in compiled C code that can be called directly from MATLAB through a MEX interface. 

The second issue that influences speed is the number of quadrature points chosen, and the considerable amount 
of bookkeeping which must be performed if the orbits are propagated and stored before evaluating the cost and 
constraints. For the examples presented in the next few sections, we used a simple trapezoid rule as opposed to a 
complex quadrature rule, in the summation of the discrete cost function. The number of points in the summation of 
the cost function was chosen to maximize the speed, while yielding an acceptable approximation to the exact, integral 
form of the cost function. 

The second issue to consider is the organization of the code. To enable convenient application to new problems, it 
is helpful to isolate code that is specific to the user-defined function. By separating and organizing the code carefully, 
it is relatively simple to provide a new function containing cost and derivatives terms. For example, for the MMS 
example in the next section, there are only 16 lines of code that are specific to the MMS cost function, and they are 
all located in a single MATLAB function that is easy to modify to solve a new problem. 

Now that we have discussed a few of the important issues in implementation, let’s take a look at a specific mission 

y • 1 • n IT 1 , >n 1 ■ / TV TCI'N 
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Application: MMS 

In the preceding portions of this paper we posed a general parameter optimization approach, and an approach to its 
implementation, to find optimal initial conditions of spacecraft formations. In this section, we’ll discuss an applica- 
tion of the method to the MMS mission. We’ll begin with a brief overview of the MMS mission. Then we’ll present 
the cost and constraint functions we choose for MMS, and derive the necessary derivatives. We conclude the section 
with a discussion of results for a particular phase of the MMS mission. 

The MMS mission is one of several missions in NASA’s Solar Terrestrial Probes (STP) program. MMS will employ 
a four spacecraft formation to make fundamental advancements in our understanding of the Earth’s magnetosphere 
and its dynamic interaction with the solar wind. MMS will not be the first mission to use multiple spacecraft to 
study magnetospheric dynamics. The Cluster II mission was successfully launched by the European Space Agency 
(ESA) in the summer of 2000 and has already provided fascinating results on magnetospheric dynamics. There are 
several phases in the MMS mission. In this work we focus on Phase I, which consists of a highly elliptic reference 
orbit with a perigee radius of 1.2 R e and an apogee radius of 12 R e . The mission design objective is to provide a near 
regular tetrahedron formation between true anomalies of 160° and 200° . There are many ways to pose a set of cost 
and constraint functions to achieve this goal and in the next section we’ll discuss the method we choose. There have 
been many important contributions to the literature on tetrahedron formations and magnetospheric missions, and a 
detailed discussion of all the past efforts is beyond the scope of this work. References [1-7,9-14] contain information 
on the science of the MMS mission, reference orbit design, launch window analysis, tetrahedron design, and maneuver 
design among other topics. 

MMS Cost and Constraint Functions 

In Phase I of the MMS mission, one of the mission design objectives is to provide a near regular tetrahedron, with 
sides of 10 km, within a region defined by ±20° in true anomaly about orbit apogee. This requires casting a set 
of cost and constraint functions that take into account the size of the tetrahedron, and its likeness to a regular 
tetrahedron. Let’s begin by investigating some useful properties of tetrahedrons. If we assume one of the spacecraft 
is the reference, then there are three relative position vectors that describe the relative geometry of the spacecraft 
and we’ll define them as s i*, S2fc> and S3* to be consistent with our previous definitions. Knowing these quantities, 
we can calculate the volume of the tetrahedron using 

V * = \ |Slfc ' (s 2 * X S 3 *) | (31) 

where V a stands for the volume of the actual tetrahedron formed by the spacecraft, as opposed to the desired volume 
which we discuss below. There are six unique sides in the formation, and the side lengths can be used to calculate 
the average side length, L *, using 


L* — - (sik + S2k + S3* + S4k + 55* + sek) (32) 

at time £*. Paschmann 10 presents several methods for evaluating how close a tetrahedron is to being a regular 
tetrahedron. One approach is to compare the volume of the actual tetrahedron V a , with the volume of a regular 
tetrahedron, V rj with side lengths equal to the average side length of the actual tetrahedron. The volume of a regular 
tetrahedron of side length L can be calculated using 


V, = 
r 12 


(33) 


Using V a and V r . 
spacecraft as 


Paschmann suggests an instantaneous volumetric performance metric, Q v , for a tetrahedron of 


. V k 

Qk _ ^a_ 

Qv ~ V k 


(34) 


where the superscripts “fc” indicate values at the k th point. Q v has the useful property: 0 < Q v < 1. However, it 
does not take into account the actual size of the tetrahedron. We propose a simple polynomial function, Q s , that 

c of a tetrahedron. T^*° ^A-nofonfo 0 . n-nH p A arp uspH to 
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change the shape of the function. A plot of Q s for a 10 km tetrahedron, with t\ — 4, £2 = 6, £3 = 18, and £4 = 25, 
is shown in Fig. 1. The important property of Q s , is that it is zero for tetrahedrons with an average side length of 
less than 4 km or greater than 25 km. Between the values of 6 km and 18 km, Q s is equal to one. 


oi(L m ) 


' 0 

(L* - h ) 2 (L* +ii- 24)7(4 - 4) 4 
< 1 

(L* - 4) 2 (L* - 24 + 4)7(4 - 4) 4 

0 


L*<£, 
ii < L* < 4 

c 2 <l* < 4 
4 < L* < 4 
L* > 4 


(35) 


By making a composite function involving both Q v and Q s we can create a performance metric that evaluates both 


Q S (L*) vs. L* 



Figure 1: Plot of Q S (L*) 


the size and shape of a tetrahedron. Let’s use the following continuous-time form 

P(t) = 1 - Qv(t)Qs(t) (36) 

P(t) will be zero for a regular tetrahedron with a side length between 6 and 18 km. For unacceptable tetrahedrons, 
P(t) is one. We can formulate a cost function by integrating P(t) over the region of interest, or 

J= [ P{t)dt- j (1 - Q v {t)Q s (t))dt 

Jto Jto 

where t Q is the time when the reference spacecraft is at a true anomaly of 160°, and tf is the time when the reference 
spacecraft is at a true anomaly of 200°. However, we can’t evaluate this integral exactly, so we have to approximate 
it. We use a simple trapezoid rule, as opposed to a complicated quadrature rule, where 


J » c£( 1 - Q k v Q k s )At k = c£( 1 - 


k= 1 


k=l 


% jSjfc • (s 2fc X S 3fc )| 

V 2 « 

/ , 3ik 
i=l 


Qsk)^tk 


(37) 


For MMS, the independent variables are the initial cartesian states of the four spacecraft in formation. In addition 
to minimizing the cost function above, we also must satisfy the nonlinear constraints that the semimajor axes of all 
spacecraft in formation are equal. We can write this as 


1 



Top ft 


(td 


(38) 


where a p is the semimajor axis of the p th spacecraft, ad is the desired semimajor axis of all spacecraft, r op is the 
magnitude of the initial position vector of the p th spacecraft, and v op is the magnitude of the initial velocity vector 
of the p th spacecraft. 
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In Eq. (37), we have a performance metric for MMS that is an explicit function of the relative geometry of the 
formation. In Eqs. (38), we have four nonlinear constraint functions in terms of the initial spacecraft states. Now all 
we require are the derivatives of the cost and constraint functions in order to apply the method to find an optimal 
solution. 


Derivatives of MMS Cost and Constraint Functions 

The last step to perform before solving for an optimal solution for MMS, is to determine the derivatives of the cost 
and constraints with respect to X G , s;/-, s^, s;*, and s ik . The k th term in the summation in the cost function, J k , 
has the form 

Jk = (1 - Q k MQl(sik)) (39) 

Notice that the cost function does not depend on the side rates, only on the side vectors and side lengths. Therefore 


and 


dJk 
ds a 


dh 

dSij 


= 0 


= 0 


The derivatives of the cost function with respect to Sik are 


dJ k _ 

dsik 

d Jk 

ds 2k 

dJk 
ds 3k ' 


( s lfc s 3 fc) Qs 


'\/2 L 3 Isijfe -S2A; x S 3fc | 

2 Sifc • S 2k X S 3 fc 
' \/2L 3 |sifc • s 2fc x s 3fc | 

c 2 Su_S2^XS3fc (s T fcS x ) Q k 

V2L 3 |sifc • s 2 k X s lfc | v lk 2k> 


where the “x” superscript indicates a skew symmetric matrix. 

Now let’s look at the derivatives of the cost function with respect to Si k . We start with 

ds ik ds lk [1 QvQs) Qv ds ik ds ik Qs 


We can show that 


where 


dQsk 

dL* 


= 


dJk _ _Q$dQ^_ jfifc • s 2 k x s 3fc | 
ds^ 6 dL* y/2 L * 4 


0 

4 (L* - w +ti- 2 e 2 )(L* - e 2 )/(-e 2 + 4) 4 
0 

4(L* - 2£ 3 + £ 4 )(L* - U)(L* - 4)/(4 - 4) 4 
0 


L* <£ t 
ii< L* < i 2 
£2 < L* < £3 
£3 < L* < £4 
L* >£4 


Finally, the nonzero derivatives of the semimajor axis constraints are 


da p _ da p dr, 
dr, 


f op 


op 


da p 

<9v. 


0 p ** op 


2 T 

= 2 ( — ] 

r op 


op 


dcip 0v op 2 a p t 

- ~ w ° p 


dv op v 0 p 


(40) 

(41) 

(42) 

(43) 

(44) 


(45) 

(46) 


(47) 


(48) 


(49) 


We now have all of the information to use an NLP code to find and optimal solutions. In the next section, we 
discuss the results for an optimal MMS formation. 


9 



Results of Optimization for MMS 

Recall that the primary formation flying goal for the phase of MMS we consider is to provide a near regular tetrahe- 
dron, with side lengths near 10 km, in the orbit region between 160° and 200° true anomaly. We also must satisfy 
the constraint that the semimajor axes for all spacecraft are the same. Figure 2 shows several plots that illustrate an 
optimal solution. The top figure shows the evolution of Q S) Q v , and Q s ■ Q v for one Keplerian orbit. The lower plot 
shows the evolution of the six individual sides and the average side length. We see that the quality factor is near 
one during the region of interest, which is bounded by the vertical dashed lines. The optimal solution is significantly 
better that the initial guess, which was calculated using an algorithm described in Ref. [9]. For the initial guess, 
the time averaged quality factor Q v - Q s was about 0.55. The NLP algorithm found an optimal solution where the 
time-averaged quality factor was about 0.94. For this solution, the run time was about 35 seconds, and required 
about 40 function evaluations. From inspection of the bottom plot in Fig. 2 we see that the average side length is 
on the order of 10 km. The orbit states for this solution are found in the Appendix. 


Figure 2: Results of MMS Orbit Optimization 


Quality Factor vs. Elapsed Time in Orbit 





Now let’s take a look at a more complicated example, the Laser Interferometer Space Antenna (LISA), which has 
a large number of nonlinear path constraints. 
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Application: LISA 

LISA is a NASA/ESA mission to detect and study gravitational waves from massive black hole systems and galactic 
binaries. Gravitational waves are ripples in space-time caused by massive objects. They are a prediction of general 
relativity, but are yet to be directly detected. The detection and understanding of gravitational waves can provide 
breakthroughs and refinements in current relativistic theory. For a more detailed discussion of the LISA mission, see 
Refs. [8, 15-17]. 

The nominal LISA formation consists of three spacecraft in heliocentric orbits trailing the Earth by about 20°, 
with inclinations near 1° with respect to the ecliptic plane. The mission design goals for LISA are challenging. The 
primary goal is to provide a formation that maintains a nearly equilateral triangle with sides near 5xl0 6 km for 
the entire life of the mission, which is currently about 8.5 years. This has to be achieved entirely through careful 
orbit design, as continuous feedback control of the orbits is not permitted because it will interfere with the science 
measurements. We also must ensure that the sides of the triangle remain within one percent of 5x1 0 6 km, and that 
the side rates never exceed 15 m/s. There is a secondary, and competing goal, that we keep the formation as close to 
Earth as possible for power reasons. Let’s look at how we can cast these goals and constraints in a manner consistent 
with the algorithm presented in previous sections. 

LISA Cost and Constraint Functions 

There are essentially two goals for LISA that can be formulated as a cost function. The first is the desire to keep the 
side lengths as close to 5xl0 6 km as possible. The second is to minimize the distance between the LISA formation 
and the Earth. We choose to pose the cost function in terms of the Earth distance, and cast the side length goals as 
constraints. The following cost function permits us to minimize the distance of the LISA formation from Earth: 


J = cY](rn +r 2 k +r 3 k) 

k-1 


(50) 


1 /2 

where = (x?fc + yf k + zf k ) 7 and so on. To ensure that the sides of the formation never exceed a one percent 
variation from the desired side length, we can apply the following set of i ■ k path constraints. 

(s ik - L) 2 < gi{t k f (51) 


where L = 5x1 0 6 and 


9\{tk) = 


-2000t fc + 50,000/cm tfc < 5 yrs 
40, 000/cm tk > 5 yrs 


(52) 


where tk is the time at the k th point, in units of years. The function pi(tfc) above, and p 2 (tfc) below, were suggested 
by Sweetser. 17 pi (tfc) and p 2 (tfc) are chosen so that the range and range rate requirements are still satisfied with 
expected orbit insertion errors. The following set of i ■ k constraints ensures that the side rates never exceed 15 m/s: 


(Sik) 2 < 92{tk) 2 


where 

92(t) = | 

Finally, we require that the ecliptic inclination, 


--tk + 15m/ s t fc < 5 yrs 
5 

12 m/s tfc > 5 yrs 

z ec , of each orbit is less than one degree. 


(53) 

(54) 


< 1 ° 


Let’s reformulate this constraint so that the derivatives are simpler. 

COS f ec = Z T h ec > COS 1° 


where 

z = (0 0 ll T 


(55) 

(56) 

(57) 
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hec = r ec x v ec (58) 

and r ec and v ec are the initial position and velocity of a LISA spacecraft, in the Sun-centered mean ecliptic of J2000 
frame, and are given by 

r ec = r fi/s + Rr ofc ; (59) 

v ec = V£ /5 + Rv o/c ; (60) 

where r E /s an d v E /s are position and velocity of the Earth with respect to the Sun, at the initial epoch, and R 
is the rotation matrix from the Earth’s mean equator of J2000 to the mean ecliptic of J2000. 

The derivatives of the cost and constraint functions are presented, without derivation, in Table 1. The leftmost 
column in the table indicates the function, and the rows are the derivatives with respect to the term that appears 
in the column labels. The superscript “x” indicates the skew symmetric matrix. For derivatives with respect to the 
initial cartesian states, the u p n subscript is a dummy variable associated with the state of the p th spacecraft. 


Table 1: Derivatives of LISA Cost and Constraint Functions 


Function 


OSih 


dv op 

nk 



r 7 ! 

r'i' 

dcJ2(rik+r 2k +r 3k ) 

0 

0 

—A p(t fc ,t 0 ) 

— B p (t k ,t 0 ) 

fc= 1 



r pk 

r pk 

d(s ik - L) 2 

2 {s ik - L) 

0 

0 

0 

d(s ik ) 2 

0 

2sifc 

0 

0 

d(x T )h ec 

0 

0 

zT (“ V E/S R - Rv °p) 

zT ( r E/S R + Rr op) 


Let’s discuss some practical considerations before moving on to the results of the LISA optimization. The current 
launch date for LISA is Jan. 2015. The initial epoch we use is 01 Jan. 2015 12:00:00.000 UTC. The force model 
used in the LISA optimization includes point masses from all planets and the Earth’s moon. The coordinate system 
is Earth’s mean J2000 equator. However, the converged solution is presented in the Appendix in the sun-centered 
mean ecliptic of J2000. Now let’s look at details of a representative initial guess and optimal solution for LISA. 

Results of Optimization for LISA 

Figure 3 contains plots illustrating characteristics of the initial guess and converged solution for an 8.5 year LISA 
trajectory. The plots associated with the initial guess are in the left column, the plots associated with the converged 
solution are in the right column. The independent variable for each plot is the elapsed mission time in years. The 
first row of plots shows the side lengths between spacecraft, or s*. The second row of plots shows the side rates, or 
si . Finally, the third row contains plots illustrating the evolution of the spacecraft-to-Earth distance. 

Let’s begin by taking a closer look at some properties of the initial guess. The first point we see is that the initial 
guess is not a feasible solution. We see that the maximum variation in range between any two LISA spacecraft is 
around 700,000 km. This is over an order of magnitude larger than the desired maximum range variation of 50,000 
km. Similarly, the maximum range rate for the initial guess is around 140 m/s, which is about an order of magnitude 
higher than the 15 m/s constraint. From inspection of the plots, we see that the initial guess is unstable. It turns 
out that the instability is due to the fact that the initial guess has been placed too close to the Earth. We see that 
the distance from Earth to the formation varies from a minimum of about 31xl0 6 km, to about 53xl0 6 km. 

The optimal solution is found to have significantly improved properties over the initial guess. Table 2 contains 
a comparison of relevant statistics for the two cases. We see that the maximum variation in range between any two 
spacecraft is 50,000 km for the converged solution. This indicates that the range constraint is active at the solution. 
The range rate between any two spacecraft is 13.7 m/s. The penalty we have paid to satisfy these constraints, is that 
the formation was moved farther away from Earth. The minimum and maximum distance from any LISA spacecraft 
to the Earth, is about 49xl0 6 km and 64xl0 6 km respectively. 

In summary, we have seen that the method we present here can solve a complex mission design problem, with 
multiple real-world constraints, in a deep-space flight regime. We began with a set of goals and requirements for 
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Figure 3: Initial Guess and Converged Solution for LISA Optimization 







the LISA mission, and cast them mathematically as a set of cost and constraint functions. Next, we presented the 
derivatives of the cost and constraint functions with respect to the formation geometry. Comparing the initial guess 
to the converged solution, indicates that we can satisfy the current mission requirements for LISA, and that we don’t 
need a feasible guess to find a solution. Now let’s look at some general conclusions we can draw from this work. 


Table 2: Statistics for LISA Initial Guess and Converged Solution 


Trajectory 

max(\sik |) (m/s) 

max(sik) (10 6 km) 

min(sik) (10^ km) 

max ( rik ) ( 10 b km) 

Initial Guess 

140 

5.7 

4.4 

53 

Converged Solution 

13.7 

5.04 

4.95 

64 
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Conclusions 


We began this work with a long list of objectives that is motivated by the need for a practical, efficient method to 
solve complex, real-world guidance problems to support a diverse set of formation flying missions. The approach we 
developed and presented meets many, if not all of our initial goals. The approach is non-linear and not restricted to 
specific flight regimes or small inter-spacecraft separations. By carefully implementing the method, we can minimize 
the amount of work that is required to solve new problems. For example, the entire MATLAB implementation of the 
cost and constraint functions for the MMS example, neglecting code for propagation and the NLP solver, consists 
of only 354 lines of code. Only 36 lines are specific to the MMS problem. The reason so few lines of code are 
specific to the mission problem is due to the fact that we assumed the cost function is an explicit function of the 
spacecraft relative positions, velocities, ranges, and range rates. This simple assumption allows us to perform half of 
the derivation of the derivatives a-priori in a general way, and implement them in software only once. This reduces 
the work that must be performed to solve new analysis problems. Furthermore, the derivatives that can be calculated 
a-priori are the derivatives of the dynamics with respect to the initial conditions. The derivatives that are problem 
dependent, are often simply derivatives of algebraic functions with respect to the formation geometry. 

We illustrated the applicability of the approach to two different missions: MMS and LISA. Both of these missions 
have a complex set of cost and constraint functions that are explicitly dependent upon the formation geometry. We 
presented optimal orbit solutions for both missions. 
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Appendix: Optimal Orbit States for MMS and LISA 


Table 3: MMS Results: State Information for a Representative, Phase I, 10 km Tetrahedron Formation 


Property 

MMS1 

MMS2 

MMS3 

MMS4 

a (km) 

42095.7 

42095.7 

42095.7 

42095.7 

e 

0.81818 

0.81798411 

0.81799342 

0.81800317 

i (deg.) 

27.8 

27.800078 

27.801283 

27.798145 

w (deg.) 

15.000001 

15.008968 

15.002286 

14.982171 

Q (deg.) 

0 

359.99962 

359.98943 

0.014249285 

v (deg.) 

180 

179.99657 

180.00303 

180.00225 


Table 4: LISA Results in Heliocentric Mean Ecliptic J2000, 01 Jan 2015 12:00:00.000 UTC 


OE 

LISA1 

LISA2 

LISA3 

LISA-Center 

a (km) 

149457958 

149457278 

149457372 

149457536 

e 

0.018265676 

0.00949038651 

0.00887263317 

0.0086422348 

i (deg.) 

0.890058699 

0.993559828 

0.993293398 

0.0663052608 

uj (deg.) 

90.6009132 

32.5150215 

148.24379 

272.489236 

(deg.) 

3.64713792 

126.757235 

240.28666 

182.897519 

v (deg.) 

343.945562 

278.418085 

51.016237 

-16.9098743 
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